1 Introduction and notes

2 Descriptive tables and figures

2.1 Measles cases and estimated vaccination coverage

Number of measles cases and estimated vaccination coverage for the 16 German states from 2005-2007. MMR1 (MMR2) is the estimated vaccination coverage for at least 1 (2) MMR vaccine, estimated from the school entry examinations. Note this partially reproduces Table 1 from Herzog et al 2011.

                                State Population Total.Cases  MMR1  MMR2
1             Baden-Wuerttemberg (BW) 10,738,753         162 90.0% 75.6%
2                        Bavaria (BY) 12,492,658         606 88.7% 73.2%
3                         Berlin (BE)  3,404,037         104 90.0% 80.2%
4                    Brandenburg (BB)  2,547,772          18 93.9% 86.9%
5                         Bremen (HB)    663,979           4 88.4% 71.9%
6                        Hamburg (HH)  1,754,182          29 90.0% 80.5%
7                          Hesse (HE)  6,075,359         336 91.2% 78.1%
8  Mecklenburg-Western Pomerania (MV)  1,693,754           4 93.6% 88.0%
9                   Lower Saxony (NI)  7,982,685         144 91.2% 78.0%
10        North Rhine-Westphalia (NW) 18,028,745       2,036 89.7% 76.9%
11          Rhineland-Palatinate (RP)  4,052,860          85 90.8% 77.3%
12                      Saarland (SL)  1,043,167           0 91.0% 81.8%
13                        Saxony (SN)  4,249,774          18 94.3% 82.4%
14                 Saxony-Anhalt (ST)  2,441,787          12 94.1% 86.5%
15            Schleswig-Holstein (SH)  2,834,254          89 89.9% 79.3%
16                     Thuringia (TH)  2,311,140           8 94.8% 85.9%

2.2 Summary figures

2.3 Maps of the estimated proportion of individuals with at least one and two MMR vaccines.

2.4 Observed bi-weekly incidence of reported cases:

3 Fit hierarchical model in Stan

3.1 Get prior

priorch <- function(x,q1,q2,p1,p2){
  (p1-pbeta(q1,x[1],x[2]))^2 + (p2-pbeta(q2,x[1],x[2]))^2 }
# p1 and p2 define the upper and lower levels of the interval
p1 <- 0.05; p2 <- 0.95 
# q1 and q2 define the bounds of the interval
q1 <- 0.6; q2 <- 0.95 
## together, we interpret these values to say:
## we want 90% (p2-p1) of the mass to be between q1 and q2

## solve and get the parameters for a beta distribution 90% of mass between 0.6 and 0.95
opt <- optim(par=c(1,1),fn=priorch,q1=q1,q2=q2,p1=p1,p2=p2)
opt$par
[1] 10.004052  2.449006

3.2 Fit model and assess converged

measlesdat=list(I=16, # number of areas
                Ts=78, # number of time points
                Nobs=16*78, # total number of observations
                sin1=sin3, cos1=cos3, # inputs for seasonal terms
                y=t(Y), # observed cases dim I x Ts
                logNi=log(pop/sum(pop)), # log pop fraction
                X=X1, # vaccination coverage
                N=sum(pop), # total pop
                a=c(10, 2.5) # parameters specifying prior on phi
                ) 


pmt=proc.time()
fit <- stan(
        model_code = mod_ENre,
        chains=4, iter=4e3, 
        control=list(adapt_delta=0.99, max_treedepth = 15),
        seed=47,
        data = measlesdat)
In file included from C:/Users/lfisher/Documents/R/win-library/3.5/BH/include/boost/config.hpp:39:0,
                 from C:/Users/lfisher/Documents/R/win-library/3.5/BH/include/boost/math/tools/config.hpp:13,
                 from C:/Users/lfisher/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/core/var.hpp:7,
                 from C:/Users/lfisher/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
                 from C:/Users/lfisher/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/core.hpp:12,
                 from C:/Users/lfisher/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/mat.hpp:4,
                 from C:/Users/lfisher/Documents/R/win-library/3.5/StanHeaders/include/stan/math.hpp:4,
                 from C:/Users/lfisher/Documents/R/win-library/3.5/StanHeaders/include/src/stan/model/model_header.hpp:4,
                 from file2cec7ae77eb9.cpp:8:
C:/Users/lfisher/Documents/R/win-library/3.5/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
 #  define BOOST_NO_CXX11_RVALUE_REFERENCES
 ^
<command-line>:0:0: note: this is the location of the previous definition
(proc.time()-pmt)/60
       user      system     elapsed 
 0.04483333  0.03866667 17.28033333 
print(fit, pars=c("alpha_ar", "phi",  'alpha_en', 'gamma', 'delta', 'sigma_ar', 'sigma_en', "lp__"), probs=c(.025,.5,.975)) # 
Inference for Stan model: c4c1ff19f358c6c9ee8fb5649733aef4.
4 chains, each with iter=4000; warmup=2000; thin=1; 
post-warmup draws per chain=2000, total post-warmup draws=8000.

            mean se_mean   sd    2.5%     50%   97.5% n_eff Rhat
alpha_ar    0.85    0.01 0.50   -0.26    0.91    1.66  1166    1
phi         0.89    0.00 0.09    0.66    0.92    0.99  2508    1
alpha_en    3.47    0.01 0.43    2.54    3.52    4.15  1956    1
gamma       0.71    0.00 0.08    0.55    0.71    0.86  8000    1
delta      -0.20    0.00 0.08   -0.36   -0.20   -0.04  8000    1
sigma_ar    0.74    0.01 0.35    0.28    0.66    1.61   685    1
sigma_en    0.55    0.00 0.17    0.28    0.52    0.95  1963    1
lp__     8516.14    0.17 5.65 8504.09 8516.47 8526.34  1153    1

Samples were drawn using NUTS(diag_e) at Mon Apr 29 16:18:41 2019.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).
print(fit, pars=c("alpha_ar", 're_ar', 'sigma_ar', "lp__"), probs=c(.025,.5,.975)) # 
Inference for Stan model: c4c1ff19f358c6c9ee8fb5649733aef4.
4 chains, each with iter=4000; warmup=2000; thin=1; 
post-warmup draws per chain=2000, total post-warmup draws=8000.

             mean se_mean   sd    2.5%     50%   97.5% n_eff Rhat
alpha_ar     0.85    0.01 0.50   -0.26    0.91    1.66  1166    1
re_ar[1]     0.43    0.01 0.31   -0.07    0.38    1.18   795    1
re_ar[2]     0.55    0.01 0.31    0.08    0.50    1.29   679    1
re_ar[3]     0.30    0.01 0.32   -0.22    0.26    1.07   755    1
re_ar[4]     0.30    0.01 0.43   -0.49    0.28    1.21  1838    1
re_ar[5]    -0.62    0.01 0.70   -2.28   -0.51    0.46  3029    1
re_ar[6]    -0.80    0.01 0.55   -2.14   -0.71    0.03  3129    1
re_ar[7]     0.64    0.01 0.30    0.19    0.60    1.38   709    1
re_ar[8]    -0.26    0.01 0.65   -1.78   -0.20    0.86  5237    1
re_ar[9]     0.20    0.01 0.32   -0.32    0.16    0.95  1017    1
re_ar[10]    0.76    0.01 0.30    0.32    0.71    1.49   671    1
re_ar[11]    0.20    0.01 0.34   -0.37    0.16    0.96  1094    1
re_ar[12]    0.02    0.01 0.87   -1.69    0.01    1.80  6090    1
re_ar[13]   -0.44    0.01 0.54   -1.74   -0.37    0.48  4084    1
re_ar[14]   -0.93    0.01 0.70   -2.60   -0.81    0.10  2232    1
re_ar[15]    0.37    0.01 0.32   -0.15    0.33    1.12   959    1
re_ar[16]   -0.74    0.01 0.71   -2.46   -0.62    0.33  2522    1
sigma_ar     0.74    0.01 0.35    0.28    0.66    1.61   685    1
lp__      8516.14    0.17 5.65 8504.09 8516.47 8526.34  1153    1

Samples were drawn using NUTS(diag_e) at Mon Apr 29 16:18:41 2019.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).
print(fit, pars=c("alpha_en", 're_en', 'sigma_en', "lp__"), probs=c(.025,.5,.975)) # 
Inference for Stan model: c4c1ff19f358c6c9ee8fb5649733aef4.
4 chains, each with iter=4000; warmup=2000; thin=1; 
post-warmup draws per chain=2000, total post-warmup draws=8000.

             mean se_mean   sd    2.5%     50%   97.5% n_eff Rhat
alpha_en     3.47    0.01 0.43    2.54    3.52    4.15  1956    1
re_en[1]    -0.26    0.00 0.22   -0.70   -0.25    0.17  3013    1
re_en[2]     0.34    0.00 0.22   -0.08    0.33    0.80  2408    1
re_en[3]     0.58    0.00 0.25    0.11    0.57    1.09  2588    1
re_en[4]    -0.32    0.00 0.33   -1.02   -0.31    0.30  8000    1
re_en[5]    -0.10    0.00 0.36   -0.84   -0.09    0.60  8000    1
re_en[6]     0.61    0.01 0.27    0.09    0.60    1.16  2505    1
re_en[7]     0.00    0.00 0.25   -0.51    0.00    0.50  3438    1
re_en[8]    -0.56    0.01 0.39   -1.40   -0.54    0.13  4267    1
re_en[9]     0.37    0.00 0.22   -0.06    0.36    0.81  2547    1
re_en[10]   -0.01    0.00 0.21   -0.42   -0.02    0.41  2514    1
re_en[11]    0.45    0.00 0.25   -0.04    0.45    0.95  2989    1
re_en[12]   -0.85    0.01 0.50   -2.01   -0.79   -0.04  3283    1
re_en[13]   -0.41    0.00 0.29   -1.01   -0.40    0.13  3819    1
re_en[14]   -0.04    0.00 0.29   -0.62   -0.03    0.52  4432    1
re_en[15]    0.50    0.01 0.26    0.00    0.49    1.02  2688    1
re_en[16]   -0.24    0.00 0.32   -0.90   -0.23    0.35  8000    1
sigma_en     0.55    0.00 0.17    0.28    0.52    0.95  1963    1
lp__      8516.14    0.17 5.65 8504.09 8516.47 8526.34  1153    1

Samples were drawn using NUTS(diag_e) at Mon Apr 29 16:18:41 2019.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).
traceplot(fit, pars=c("alpha_ar", "phi",  'alpha_en', 'gamma', 'delta','sigma_ar', 'sigma_en', "lp__")) #

pairs(fit, pars=c("alpha_ar", "phi",  'alpha_en', 'gamma', 'delta', 'sigma_ar', 'sigma_en', "lp__")) #

3.3 Parameter estimates

            50%   2.5%  97.5%
alpha_ar  0.911 -0.256  1.656
phi       0.917  0.659  0.989
alpha_en  3.525  2.543  4.155
gamma     0.708  0.551  0.860
delta    -0.200 -0.356 -0.043
sigma_ar  0.663  0.277  1.607
sigma_en  0.524  0.279  0.955
R0        2.488  0.774  5.237

3.4 Examine model fits

3.5 Estimates of random effects

3.6 Examine posterior distributions

3.7 Plot Pearson residuals

4 Additional analyses

4.1 Comparison of MLE vs Bayesian analysis

4.1.1 MLE estimates

loglik_eco=function(param, y, xi, n){
        v=expit(param[2])
        l=0
        nwks=dim(y)[2]
        ni=dim(y)[1]
        for(i in 1:ni){
                for(t in 2:nwks){
                        lam=(exp(param[1])*y[i, t-1] + (n[i]/sum(n))*exp(param[3] + param[4]*sin(2*pi*t/26) + param[5]*cos(2*pi*t/26))  )
                        l=l-dpois(y[i, t], lambda=(1-v*xi[i])*lam, log=T)
                }
        }
        return(l)
}
  parameter         Est         SE
1  alpha_ar  2.02928493 0.07239919
2       phi  0.98867494 0.68640739
3  alpha_en  4.06273036 0.10100270
4     gamma  0.69281062 0.08418688
5     delta -0.07449201 0.08120128

4.1.2 Estimates from fixed effects Stan model

measlesdat2=measlesdat
measlesdat2$a=c(1, 1) # put a flat prior on phi

pmt=proc.time()
fit_fe <- stan(
        model_code = mod_fe,
        chains=4, iter=4e3, # warmup = 1e3, # save_dso = T,  # warmup=7e3, # refresh=F, #  #
        control=list(adapt_delta=0.99, max_treedepth = 15),
        data = measlesdat2)
In file included from C:/Users/lfisher/Documents/R/win-library/3.5/BH/include/boost/config.hpp:39:0,
                 from C:/Users/lfisher/Documents/R/win-library/3.5/BH/include/boost/math/tools/config.hpp:13,
                 from C:/Users/lfisher/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/core/var.hpp:7,
                 from C:/Users/lfisher/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
                 from C:/Users/lfisher/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/core.hpp:12,
                 from C:/Users/lfisher/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/mat.hpp:4,
                 from C:/Users/lfisher/Documents/R/win-library/3.5/StanHeaders/include/stan/math.hpp:4,
                 from C:/Users/lfisher/Documents/R/win-library/3.5/StanHeaders/include/src/stan/model/model_header.hpp:4,
                 from file29ac349f2e0d.cpp:8:
C:/Users/lfisher/Documents/R/win-library/3.5/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
 #  define BOOST_NO_CXX11_RVALUE_REFERENCES
 ^
<command-line>:0:0: note: this is the location of the previous definition
(proc.time()-pmt)/60
      user     system    elapsed 
0.06216667 0.04533333 3.48383333 
print(fit_fe, pars=c("alpha_ar", "phi",  'alpha_en', 'gamma', 'delta'), probs=c(.025,.5,.975), digits_summary = 3) # 
Inference for Stan model: 8a16ff8514700ac19fc4f3422c38a24a.
4 chains, each with iter=4000; warmup=2000; thin=1; 
post-warmup draws per chain=2000, total post-warmup draws=8000.

           mean se_mean    sd   2.5%    50% 97.5% n_eff  Rhat
alpha_ar  2.066   0.002 0.057  1.918  2.080 2.141  1380 1.003
phi       0.993   0.000 0.007  0.975  0.996 1.000  1307 1.003
alpha_en  4.120   0.002 0.087  3.922  4.126 4.273  2068 1.002
gamma     0.698   0.001 0.084  0.534  0.699 0.863  5049 1.001
delta    -0.074   0.001 0.080 -0.232 -0.074 0.082  4526 1.000

Samples were drawn using NUTS(diag_e) at Thu May 02 14:29:44 2019.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

4.2 Sensitivity analysis - measles analysis with non-informative priors

measlesdat=list(I=16, # number of areas
                Ts=78, # number of time points
                Nobs=16*78, # total number of observations
                sin1=sin3, cos1=cos3, # inputs for seasonal terms
                y=t(Y), # observed cases dim I x Ts
                logNi=log(pop/sum(pop)), # log pop fraction
                X=X1, # vaccination coverage
                N=sum(pop), # total pop
                a=c(1, 1) # parameters specifying prior on phi
                ) 


pmt=proc.time()
fit_flat <- stan(
        model_code = mod_ENre,
        chains=4, iter=4e3, # warmup = 1e3, # save_dso = T,  # warmup=7e3, # refresh=F, #  #
        control=list(adapt_delta=0.99, max_treedepth = 15),
        data = measlesdat)
(proc.time()-pmt)/60
       user      system     elapsed 
 0.01650000  0.02766667 21.40100000 
print(fit_flat, pars=c("alpha_ar", "phi",  'alpha_en', 'gamma', 'delta', 'sigma_ar', 'sigma_en', "lp__"), probs=c(.025,.5,.975)) # 
Inference for Stan model: c4c1ff19f358c6c9ee8fb5649733aef4.
4 chains, each with iter=4000; warmup=2000; thin=1; 
post-warmup draws per chain=2000, total post-warmup draws=8000.

            mean se_mean   sd    2.5%     50%   97.5% n_eff Rhat
alpha_ar    0.97    0.03 0.77   -0.91    1.24    1.85   687    1
phi         0.87    0.01 0.22    0.15    0.97    1.00   892    1
alpha_en    3.58    0.02 0.70    1.83    3.85    4.32   871    1
gamma       0.71    0.00 0.08    0.55    0.71    0.87  7386    1
delta      -0.20    0.00 0.08   -0.36   -0.20   -0.03  6306    1
sigma_ar    0.69    0.01 0.33    0.24    0.62    1.53   772    1
sigma_en    0.51    0.00 0.18    0.24    0.49    0.95  1567    1
lp__     8521.75    0.21 6.10 8509.33 8522.09 8532.99   841    1

Samples were drawn using NUTS(diag_e) at Mon Apr 29 16:44:34 2019.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).
print(fit_flat, pars=c("alpha_ar", 're_ar', 'sigma_ar', "lp__"), probs=c(.025,.5,.975)) # 
Inference for Stan model: c4c1ff19f358c6c9ee8fb5649733aef4.
4 chains, each with iter=4000; warmup=2000; thin=1; 
post-warmup draws per chain=2000, total post-warmup draws=8000.

             mean se_mean   sd    2.5%     50%   97.5% n_eff Rhat
alpha_ar     0.97    0.03 0.77   -0.91    1.24    1.85   687    1
re_ar[1]     0.38    0.01 0.31   -0.11    0.34    1.11   792    1
re_ar[2]     0.50    0.01 0.31    0.03    0.45    1.22   707    1
re_ar[3]     0.27    0.01 0.32   -0.24    0.23    1.01   849    1
re_ar[4]     0.30    0.01 0.42   -0.51    0.28    1.15  1992    1
re_ar[5]    -0.58    0.01 0.66   -2.25   -0.47    0.45  3012    1
re_ar[6]    -0.75    0.01 0.51   -1.91   -0.68    0.03  3393    1
re_ar[7]     0.61    0.01 0.29    0.17    0.57    1.28   793    1
re_ar[8]    -0.23    0.01 0.61   -1.62   -0.16    0.80  4131    1
re_ar[9]     0.18    0.01 0.31   -0.33    0.14    0.88   959    1
re_ar[10]    0.72    0.01 0.29    0.28    0.67    1.41   728    1
re_ar[11]    0.17    0.01 0.32   -0.38    0.14    0.90  1108    1
re_ar[12]   -0.01    0.01 0.75   -1.56   -0.01    1.50  6256    1
re_ar[13]   -0.39    0.01 0.52   -1.59   -0.32    0.47  4406    1
re_ar[14]   -0.84    0.02 0.68   -2.47   -0.73    0.15  1849    1
re_ar[15]    0.33    0.01 0.32   -0.19    0.29    1.07   894    1
re_ar[16]   -0.66    0.01 0.66   -2.24   -0.55    0.32  1964    1
sigma_ar     0.69    0.01 0.33    0.24    0.62    1.53   772    1
lp__      8521.75    0.21 6.10 8509.33 8522.09 8532.99   841    1

Samples were drawn using NUTS(diag_e) at Mon Apr 29 16:44:34 2019.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).
print(fit_flat, pars=c("alpha_en", 're_en', 'sigma_en', "lp__"), probs=c(.025,.5,.975)) # 
Inference for Stan model: c4c1ff19f358c6c9ee8fb5649733aef4.
4 chains, each with iter=4000; warmup=2000; thin=1; 
post-warmup draws per chain=2000, total post-warmup draws=8000.

             mean se_mean   sd    2.5%     50%   97.5% n_eff Rhat
alpha_en     3.58    0.02 0.70    1.83    3.85    4.32   871    1
re_en[1]    -0.27    0.00 0.23   -0.72   -0.27    0.19  2620    1
re_en[2]     0.30    0.01 0.23   -0.13    0.28    0.78  1654    1
re_en[3]     0.54    0.01 0.26    0.07    0.53    1.05  2166    1
re_en[4]    -0.27    0.01 0.33   -0.97   -0.25    0.33  3655    1
re_en[5]    -0.11    0.00 0.36   -0.86   -0.10    0.58  5456    1
re_en[6]     0.56    0.01 0.28    0.02    0.55    1.12  2166    1
re_en[7]     0.00    0.00 0.25   -0.49    0.00    0.48  3239    1
re_en[8]    -0.51    0.01 0.39   -1.39   -0.47    0.15  3284    1
re_en[9]     0.35    0.00 0.22   -0.07    0.35    0.80  2379    1
re_en[10]   -0.03    0.00 0.21   -0.42   -0.04    0.39  2050    1
re_en[11]    0.42    0.01 0.26   -0.06    0.41    0.94  2378    1
re_en[12]   -0.79    0.01 0.50   -1.97   -0.72   -0.01  3016    1
re_en[13]   -0.36    0.01 0.30   -0.98   -0.34    0.18  3299    1
re_en[14]   -0.01    0.01 0.29   -0.60    0.00    0.53  3066    1
re_en[15]    0.46    0.01 0.26   -0.03    0.45    1.01  2159    1
re_en[16]   -0.20    0.01 0.32   -0.88   -0.18    0.40  3619    1
sigma_en     0.51    0.00 0.18    0.24    0.49    0.95  1567    1
lp__      8521.75    0.21 6.10 8509.33 8522.09 8532.99   841    1

Samples were drawn using NUTS(diag_e) at Mon Apr 29 16:44:34 2019.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).
traceplot(fit_flat, pars=c("alpha_ar", "phi",  'alpha_en', 'gamma', 'delta','sigma_ar', 'sigma_en', "lp__")) #

pairs(fit_flat, pars=c("alpha_ar", "phi",  'alpha_en', 'gamma', 'delta', 'sigma_ar', 'sigma_en', "lp__")) #

4.2.1 Examine model fits

4.2.2 Estimates of random effects

4.2.3 Examine posterior distributions

5 Stan Models

Below, we print the code to define the Stan models fit in this document

mod_ENre='
data{
int<lower=1> I;          // number of areas
int<lower=1> Ts;         // number of times
int<lower=1> Nobs;    // number of total observations
int y[I, Ts] ;  // IxTs matrix of observations
real logNi[I];          // log population proportions 
vector[I] X;            // vector of Xi vacc cover
vector[Ts] sin1;           // vector of sin
vector[Ts] cos1;           // vector of cos
vector[2] a; //vector of parameters for the beta prior on phi
}
parameters{
real alpha_ar;  // AR component - intercept
real alpha_en;  // EN component - intercept
real gamma;     // sin component
real delta;     // cosine component 
real<lower=0, upper=1> phi;  // vaccine effect
vector[I] re_en; // EN random effects
vector[I] re_ar; // AR random effects
real<lower=0> sigma_ar; 
real<lower=0> sigma_en; 
}
model{
for(i in 1:I){
re_ar[i]~normal(0, sigma_ar);
re_en[i]~normal(0, sigma_en);
for(t in 2:Ts){
real muit;
muit = exp(alpha_ar+re_ar[i])*y[i, t-1] + exp(logNi[i] + alpha_en + re_en[i] + gamma*sin1[t] + delta*cos1[t]);
y[i, t] ~ poisson((1-phi*X[i])*muit);
}
}

//priors
alpha_ar ~ normal(0, 5);
alpha_en ~ normal(0, 5);
gamma ~ normal(0, 10);
delta ~ normal(0, 10);
phi ~ beta(a[1], a[2]);

}
generated quantities{
vector[Nobs-I] log_lik;
vector[I] ar;
vector[I] en;
vector[I] r; 
for(i in 1:I){
r[i] = exp(alpha_ar + re_ar[i])*(1-phi*X[i]); 
ar[i] = exp(alpha_ar + re_ar[i]);
en[i] = exp(alpha_en + re_en[i]);
for(t in 2:Ts){
real muit;
muit=exp(alpha_ar + re_ar[i])*y[i, t-1] + exp(logNi[i] + alpha_en + re_en[i] + gamma*sin1[t] + delta*cos1[t]);
log_lik[(i-1)*(Ts-1)+t-1 ] = poisson_lpmf(y[i, t] | (1-phi*X[i])*muit);
}
}
}
'

### fixed effects model - to the MLE 
mod_fe='
data{
int<lower=1> I;          // number of areas
int<lower=1> Ts;         // number of times
int<lower=1> Nobs;    // number of total observations
int y[I, Ts] ;  // IxTs matrix of observations
real logNi[I];          // log population sizes
vector[I] X;            // vector of Xi vacc cover
vector[Ts] sin1;           // vector of sin
vector[Ts] cos1;           // vector of cos
vector[2] a; //vector of parameters for the beta prior on phi
}
parameters{
real alpha_ar;  // AR component - intercept
real alpha_en;  // EN component - intercept
real gamma;     // sin component
real delta;     // cosine component 
real<lower=0, upper=1> phi;  // vaccine effect
}
model{
for(i in 1:I){
for(t in 2:Ts){
real muit;
muit = exp(alpha_ar)*y[i, t-1] + exp(logNi[i] + alpha_en + gamma*sin1[t] + delta*cos1[t]);
y[i, t] ~ poisson((1-phi*X[i])*muit);
}
}

//priors
alpha_ar ~ normal(0, 5);
alpha_en ~ normal(0, 5);
gamma ~ normal(0, 10);
delta ~ normal(0, 10);
phi ~ beta(a[1], a[2]);

}
generated quantities{
vector[Nobs-I] log_lik;
vector[I] ri; 
for(i in 1:I){
ri[i] = exp(alpha_ar)*(1-phi*X[i]); 
for(t in 2:Ts){
real muit;
muit=exp(alpha_ar)*y[i, t-1] + exp(logNi[i] + alpha_en + gamma*sin1[t] + delta*cos1[t]);
log_lik[(i-1)*(Ts-1)+t-1 ] = poisson_lpmf(y[i, t] | (1-phi*X[i])*muit);
}
}
}
'